You’ll need to download Stan to your machine for these analyses. Note: this document will save brms model fits into your working directory.
First, let’s load all the packages we’ll need.
library(bayesplot)
library(brms)
library(cowplot)
library(ggplot2)
library(grid)
library(tidyverse)
Now, read in the data from Study 2.
(
d <-
read.csv('data/exp2data.csv', header = TRUE) %>%
as_tibble()
)
Each row is a different participant in Study 2 (there were 80). Columns 1-68 are variables that summarise the whole session. Columns 69-493 refer to events in each round of the Visible condition. Columns 494-918 refer to the Hidden condition. This experiment followed a within-subjects design (i.e. participants completed both the Hidden and Visible conditions). In the counterbalancing column, 0 means that the participant played the Visible game before the Hidden game, and 1 means that the participant played the Hidden game before the Visible game.
Trim the dataset and get it in long-format.
(
d.trim <-
d %>%
select(ID, ends_with('.rounds_survived'), ends_with('.overall_num_shocks')) %>%
gather(key, value, -ID) %>%
extract(key, c("Condition", "variable"),
"(hidden|visible).(rounds_survived|overall_num_shocks)") %>%
spread(variable, value) %>%
mutate(Condition = ifelse(Condition == 'hidden', 1, 0))
)
As in Study 1, we first test whether the probability of a shock occurring differs between the two conditions. Although this experiment followed a within-subjects design, random effects are not necessary here because the occurrence of shocks are independent between the two conditions.
All analyses in this document will use the brm() function. Learn more about the brms package here. Also, see here for an explanation of why I use the 0 + intercept syntax throughout this document.
b2.1 <-
brm(data = d.trim, family = binomial,
overall_num_shocks | trials(rounds_survived) ~ 0 + intercept + Condition,
prior = c(prior(normal(0, 1), class = b)),
iter = 2000, warmup = 1000, chains = 4, cores = 4,
seed = 2113)
save(b2.1, file = "models/exp2_brmsfit1.rda")
We set the seed to a random number, to make the results reproducible. Here are the priors that were set for this model.
prior_summary(b2.1)
Let’s look at the results.
print(b2.1)
## Family: binomial
## Links: mu = logit
## Formula: overall_num_shocks | trials(rounds_survived) ~ 0 + intercept + Condition
## Data: d.trim (Number of observations: 160)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept -1.36 0.06 -1.47 -1.25 1.00 1223 1649
## Condition -0.03 0.09 -0.21 0.14 1.00 1244 1612
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Plotting the parameters.
stanplot(b2.1)
The Condition parameter is -0.03, and its 95% credible intervals cross 0, implying that the probability of a shock is no different across the two conditions. On the probability scale:
post <- posterior_samples(b2.1)
visible_prob <- inv_logit_scaled(post$b_intercept)
hidden_prob <- inv_logit_scaled(post$b_intercept + post$b_Condition)
difference <- hidden_prob - visible_prob
quantile(difference,c(0.025,0.5,0.975)) %>% round(2)
## 2.5% 50% 97.5%
## -0.03 -0.01 0.02
Trim the dataset again.
(
d.trim <-
d %>%
select(hidden.total_cattle_lost, visible.total_cattle_lost) %>%
gather(Condition, total_cattle_lost) %>%
mutate(Condition = ifelse(Condition == 'hidden.total_cattle_lost', 1, 0))
)
We now fit a model to determine if the total amount of cattle lost due to shocks varies between conditions. Again, no random effects are needed because the outcome is independent across conditions (generated stochastically by the game), despite the within-subject design.
b2.2 <-
brm(data = d.trim, family = gaussian,
total_cattle_lost ~ 0 + intercept + Condition,
prior = c(prior(normal(0, 100), class = b, coef = 'intercept'),
prior(normal(0, 5), class = b, coef = 'Condition')),
iter = 2000, warmup = 1000, chains = 4, cores = 4,
seed = 2113)
save(b2.2, file = "models/exp2_brmsfit2.rda")
The priors we used.
prior_summary(b2.2)
The results.
print(b2.2)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: total_cattle_lost ~ 0 + intercept + Condition
## Data: d.trim (Number of observations: 160)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept 51.27 2.11 47.26 55.53 1.00 2638 2655
## Condition -6.70 2.70 -11.98 -1.30 1.00 2736 2727
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 20.41 1.16 18.30 22.78 1.00 2617 2448
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Plot the parameters.
stanplot(b2.2)
On average, 7 more cattle were lost in the Hidden condition than the Visible condition.
Get a long-format data frame with binary request decisions over all rounds, for both conditions. If request == NA, player has died and been removed from the game, so we drop those rows.
(
d.trim <-
d %>%
select(ID, Group, ends_with('.player.request')) %>%
gather(key, request, -ID, -Group) %>%
extract(key, c("Condition", "round_number"),
"SurvivalGame_(Hidden|Visible)[.](.|..)[.]player[.]request") %>%
mutate(Condition = ifelse(Condition == 'Hidden', 1, 0),
round_number = as.integer(round_number)) %>%
drop_na()
)
This leaves us with 2834 request decisions.
We now fit a varying intercept and slope model, with participants nested within groups. We allow the slopes for both round number and condition to vary, to fit with the experiment’s within-subjects design (participants completed multiple rounds, in both conditions).
b2.3 <-
brm(data = d.trim, family = bernoulli,
request ~ 0 + intercept + round_number + Condition
+ (0 + intercept + round_number + Condition | Group/ID),
prior = c(prior(normal(0, 1), class = b),
prior(student_t(3, 0, 10), class = sd),
prior(lkj(1), class = cor)),
iter = 2000, warmup = 1000, chains = 4, cores = 4,
sample_prior = TRUE,
control = list(adapt_delta = 0.99),
seed = 2113)
save(b2.3, file = "models/exp2_brmsfit3.rda")
Here are the priors for the model we just fitted.
prior_summary(b2.3)
Now let’s see the results.
print(b2.3)
## Family: bernoulli
## Links: mu = logit
## Formula: request ~ 0 + intercept + round_number + Condition + (0 + intercept + round_number + Condition | Group/ID)
## Data: d.trim (Number of observations: 2834)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~Group (Number of levels: 40)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept) 0.65 0.29 0.06 1.17 1.02 288 278
## sd(round_number) 0.05 0.02 0.01 0.08 1.01 387 502
## sd(Condition) 1.65 0.31 1.07 2.30 1.00 939 1066
## cor(intercept,round_number) -0.40 0.38 -0.89 0.56 1.01 441 918
## cor(intercept,Condition) -0.25 0.33 -0.79 0.53 1.01 330 265
## cor(round_number,Condition) 0.34 0.30 -0.35 0.85 1.01 378 476
##
## ~Group:ID (Number of levels: 80)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept) 0.68 0.26 0.16 1.19 1.02 293 509
## sd(round_number) 0.04 0.02 0.00 0.08 1.02 251 779
## sd(Condition) 1.01 0.32 0.43 1.69 1.01 432 686
## cor(intercept,round_number) -0.41 0.44 -0.92 0.70 1.00 498 1404
## cor(intercept,Condition) -0.32 0.34 -0.83 0.49 1.01 509 630
## cor(round_number,Condition) 0.11 0.44 -0.78 0.84 1.01 315 602
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept -1.01 0.18 -1.34 -0.67 1.00 2916 2735
## round_number -0.04 0.01 -0.07 -0.02 1.00 2271 2736
## Condition 0.62 0.30 0.02 1.21 1.00 2793 3065
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Plotting the parameters.
stanplot(b2.3)
Do those Rhat values look okay? Larger than 1.1 is deemed an issue.
rhat(b2.3) %>%
mcmc_rhat()
How about nEff ratios? Above 0.1 is desired.
neff_ratio(b2.3) %>%
mcmc_neff()
Only a few parameters are below the threshold.
Finally, let’s see the trace plots.
plot(b2.3, ask = F)
Looks like Stan sampled efficiently.
In this model, the fixed effect of condition is 0.62, with 95% credible intervals above 0. This implies that the participants are more likely to request from their partner in the Hidden condition. Converting to the probability scale:
post <- posterior_samples(b2.3)
visible_prob <- inv_logit_scaled(post$b_intercept)
visible_prob %>%
median() %>%
round(2)
## [1] 0.27
hidden_prob <- inv_logit_scaled(post$b_intercept + post$b_Condition)
hidden_prob %>%
median() %>%
round(2)
## [1] 0.41
difference <- hidden_prob - visible_prob
quantile(difference,c(0.025,0.5,0.975)) %>% round(2)
## 2.5% 50% 97.5%
## 0.00 0.14 0.28
The absolute probability difference between the conditions is +0.14 (median), with 95% CIs above 0. Participants were more likely to request from their partner in the Hidden condition.
We compute a Bayes factor for this difference between probabilities. This is the inverse of the alternative hypothesis that the two conditions are equal.
(1 / hypothesis(b2.3, "inv_logit_scaled(intercept) - inv_logit_scaled(intercept + Condition) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 2.38
This Bayes factor implies weak support for the hypothesis that the probabilities differ across conditions.
How much variance in the outcome does this model explain? We use the bayes_R2() function to get a Bayesian estimate of R-squared.
bayes_R2(b2.3) %>% round(2)
## Estimate Est.Error Q2.5 Q97.5
## R2 0.27 0.01 0.24 0.29
This model explains around 27% of the variance.
Get a long-format data frame with binary request decisions over all rounds, for both conditions. If request == NA, player has died and been removed from the game, so we drop those rows. However, we also filter out rows in which the player was below the minimum survival threshold (64 cattle).
(
d.trim <-
d %>%
select(ID, Group, Counterbalancing, ends_with('.player.herd_size_after_shock'),
ends_with('.player.request')) %>%
gather(key, value, -ID, -Group, -Counterbalancing) %>%
extract(key, c("Condition", "round_number", "variable"),
"SurvivalGame_(Hidden|Visible)[.](.|..)[.]player[.](request|herd_size_after_shock)") %>%
spread(variable, value) %>%
mutate(Condition = ifelse(Condition == 'Hidden', 1, 0),
round_number = as.integer(round_number)) %>%
arrange(round_number, ID, Group) %>%
drop_na() %>%
filter(herd_size_after_shock >= 64)
)
This leaves us with 2368 request decisions.
We now fit a varying intercept and slope model, with participants nested within groups. Again, we allow the slopes for both round number and condition to vary.
b2.4 <-
brm(data = d.trim, family = bernoulli,
request ~ 0 + intercept + round_number + Condition
+ (0 + intercept + round_number + Condition | Group/ID),
prior = c(prior(normal(0, 1), class = b),
prior(student_t(3, 0, 10), class = sd),
prior(lkj(1), class = cor)),
iter = 2000, warmup = 1000, chains = 4, cores = 4,
sample_prior = TRUE,
control = list(adapt_delta = 0.99),
seed = 2113)
b2.4 <- add_criterion(b2.4, "loo")
save(b2.4, file = "models/exp2_brmsfit4.rda")
The priors.
prior_summary(b2.4)
The results.
print(b2.4)
## Family: bernoulli
## Links: mu = logit
## Formula: request ~ 0 + intercept + round_number + Condition + (0 + intercept + round_number + Condition | Group/ID)
## Data: d.trim (Number of observations: 2368)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~Group (Number of levels: 40)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept) 1.12 0.49 0.12 2.01 1.00 517 640
## sd(round_number) 0.08 0.03 0.01 0.15 1.01 397 491
## sd(Condition) 1.93 0.76 0.31 3.34 1.01 637 677
## cor(intercept,round_number) -0.15 0.42 -0.79 0.77 1.01 597 1505
## cor(intercept,Condition) -0.08 0.41 -0.80 0.79 1.01 520 922
## cor(round_number,Condition) 0.33 0.38 -0.60 0.91 1.00 695 952
##
## ~Group:ID (Number of levels: 80)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept) 1.53 0.36 0.85 2.25 1.00 848 1412
## sd(round_number) 0.09 0.03 0.04 0.14 1.00 600 925
## sd(Condition) 2.30 0.66 1.25 3.79 1.00 759 1842
## cor(intercept,round_number) -0.61 0.23 -0.92 -0.02 1.00 1483 2016
## cor(intercept,Condition) -0.25 0.25 -0.68 0.29 1.00 867 1497
## cor(round_number,Condition) 0.50 0.28 -0.14 0.93 1.01 506 1075
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept -1.96 0.33 -2.63 -1.33 1.00 3063 2850
## round_number -0.10 0.03 -0.16 -0.05 1.00 2159 2447
## Condition 0.18 0.49 -0.82 1.10 1.00 3298 2919
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Creating a forest plot of parameters.
stanplot(b2.4)
Viusalise Rhat values.
rhat(b2.4) %>%
mcmc_rhat()
Viusalise nEff ratios.
neff_ratio(b2.4) %>%
mcmc_neff()
Trace plots.
plot(b2.4, ask = F)
Looks like Stan sampled efficiently.
In this model, the fixed effect of condition is 0.18, with 95% credible intervals crossing 0. This implies that, at least when above the minimum survival threshold, the average probability of requesting did not differ between conditions.
Converting the fixed effects to the probability scale:
post <- posterior_samples(b2.4)
visible_prob <- inv_logit_scaled(post$b_intercept)
visible_prob %>%
median() %>%
round(2)
## [1] 0.12
hidden_prob <- inv_logit_scaled(post$b_intercept + post$b_Condition)
hidden_prob %>%
median() %>%
round(2)
## [1] 0.15
difference <- hidden_prob - visible_prob
quantile(difference,c(0.025,0.5,0.975)) %>% round(2)
## 2.5% 50% 97.5%
## -0.07 0.02 0.16
The absolute probability difference between the conditions is +0.02 (median), with 95% CIs above 0. Participants were no more likely to request from their partner in the Hidden condition.
We compute a Bayes factor for this difference between probabilities.
(1 / hypothesis(b2.4, "inv_logit_scaled(intercept) - inv_logit_scaled(intercept + Condition) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 0.32
This Bayes factor implies moderate support for the hypothesis that the probabilities are equal in each condition.
How much variance in the outcome does this model explain?
bayes_R2(b2.4) %>% round(2)
## Estimate Est.Error Q2.5 Q97.5
## R2 0.43 0.02 0.39 0.46
This model explains around 43% of the variance.
We found no effect of condition in the previous model. Is this because of order effects?
b2.4b <-
brm(data = d.trim, family = bernoulli,
request ~ 0 + intercept + round_number + Condition*Counterbalancing
+ (0 + intercept + round_number + Condition | Group/ID),
prior = c(prior(normal(0, 1), class = b),
prior(student_t(3, 0, 10), class = sd),
prior(lkj(1), class = cor)),
iter = 2000, warmup = 1000, chains = 4, cores = 4,
sample_prior = TRUE,
control = list(adapt_delta = 0.99),
seed = 2113)
save(b2.4b, file = "models/exp2_brmsfit4b.rda")
The priors.
prior_summary(b2.4b)
The results.
print(b2.4b)
## Family: bernoulli
## Links: mu = logit
## Formula: request ~ 0 + intercept + round_number + Condition * Counterbalancing + (0 + intercept + round_number + Condition | Group/ID)
## Data: d.trim (Number of observations: 2368)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~Group (Number of levels: 40)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept) 1.14 0.49 0.15 2.05 1.00 463 603
## sd(round_number) 0.08 0.03 0.01 0.14 1.01 380 583
## sd(Condition) 1.34 0.71 0.08 2.77 1.01 345 1108
## cor(intercept,round_number) -0.14 0.42 -0.78 0.78 1.00 591 1056
## cor(intercept,Condition) -0.19 0.44 -0.88 0.71 1.01 872 1968
## cor(round_number,Condition) 0.31 0.42 -0.68 0.92 1.01 788 1259
##
## ~Group:ID (Number of levels: 80)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept) 1.53 0.38 0.83 2.28 1.00 634 1346
## sd(round_number) 0.08 0.03 0.03 0.14 1.01 467 569
## sd(Condition) 2.14 0.55 1.18 3.29 1.01 693 1295
## cor(intercept,round_number) -0.60 0.25 -0.92 0.02 1.00 918 1014
## cor(intercept,Condition) -0.26 0.26 -0.69 0.33 1.01 486 935
## cor(round_number,Condition) 0.52 0.28 -0.11 0.93 1.01 399 859
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept -2.09 0.40 -2.89 -1.31 1.00 1801 2854
## round_number -0.10 0.03 -0.15 -0.05 1.00 2276 2460
## Condition -0.38 0.50 -1.41 0.56 1.00 2581 2981
## Counterbalancing 0.19 0.48 -0.79 1.10 1.00 1861 2203
## Condition:Counterbalancing 1.56 0.62 0.30 2.73 1.00 2777 2751
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Creating a forest plot of parameters.
stanplot(b2.4b)
Trace plots.
plot(b2.4b, ask = F)
Looks like Stan sampled efficiently.
In this model, the interaction parameter is 1.56, with 95% CIs above zero. This implies that there is an interaction effect. Condition only has an effect on cheating when participants played the Hidden game first.
Converting the fixed effects to the probability scale:
post <- posterior_samples(b2.4b)
visible_VisibleFirst_prob <- inv_logit_scaled(post$b_intercept)
visible_VisibleFirst_prob %>%
median() %>%
round(2)
## [1] 0.11
hidden_VisibleFirst_prob <- inv_logit_scaled(post$b_intercept + post$b_Condition)
hidden_VisibleFirst_prob %>%
median() %>%
round(2)
## [1] 0.08
difference_VisibleFirst <- hidden_VisibleFirst_prob - visible_VisibleFirst_prob
quantile(difference_VisibleFirst,c(0.025,0.5,0.975)) %>% round(2)
## 2.5% 50% 97.5%
## -0.11 -0.03 0.06
When participants play the Visible game first, the absolute probability difference between the conditions is -0.03 (median), with 95% CIs crossing 0. Participants were no more likely to request from their partner in the Hidden condition.
We compute a Bayes factor for this difference between probabilities.
(1 / hypothesis(b2.4b, "inv_logit_scaled(intercept) - inv_logit_scaled(intercept + Condition) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 0.31
This Bayes factor implies moderate support for the hypothesis that the probabilities are equal in each condition.
What about for when participants play the Hidden game first?
visible_HiddenFirst_prob <- inv_logit_scaled(post$b_intercept + post$b_Counterbalancing)
visible_HiddenFirst_prob %>%
median() %>%
round(2)
## [1] 0.13
hidden_HiddenFirst_prob <- inv_logit_scaled(post$b_intercept + post$b_Counterbalancing +
post$b_Condition + post$`b_Condition:Counterbalancing`)
hidden_HiddenFirst_prob %>%
median() %>%
round(2)
## [1] 0.34
difference_HiddenFirst <- hidden_HiddenFirst_prob - visible_HiddenFirst_prob
quantile(difference_HiddenFirst,c(0.025,0.5,0.975)) %>% round(2)
## 2.5% 50% 97.5%
## 0.01 0.20 0.43
When participants play the Hidden game first, the absolute probability difference between the conditions is +0.20 (median), with 95% CIs above 0. Participants are more likely to request from their partner when beneath the threshold in the Hidden condition.
We compute a Bayes factor for this difference between probabilities.
(1 / hypothesis(b2.4b, "inv_logit_scaled(intercept + Counterbalancing) - inv_logit_scaled(intercept + Counterbalancing + Condition + Condition:Counterbalancing) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 3.05
This Bayes factor implies moderate support for the hypothesis that the probabilities differ across conditions.
How much variance in the outcome does this model explain?
bayes_R2(b2.4b) %>% round(2)
## Estimate Est.Error Q2.5 Q97.5
## R2 0.42 0.02 0.39 0.45
This model explains around 42% of the variance.
Get long-format data frame with ‘received’ variable (i.e. how much a player received on any given round). We swap this around so it reflects how much the player gave to their partner (i.e. how much their partner received). We drop rows with NAs, since partners did not request help in that particular round. We then code whether the player gave nothing in response to the request (1) or gave at least one cattle (0).
(
d.trim <-
d %>%
select(ID, Group, ends_with('.player.received')) %>%
# swap every other value to get what the player GAVE, rather than what they RECEIVED
# swapping works because each player is next to their partner in the data frame
mutate_at(vars(ends_with(".player.received")), function(x) x[1:nrow(d) + c(1,-1)]) %>%
gather(key, received, -ID, -Group) %>%
extract(key, c("Condition", "round_number"),
"SurvivalGame_(Hidden|Visible)[.](.|..)[.]player[.]received") %>%
rename(responded = received) %>%
mutate(Condition = ifelse(Condition == 'Hidden', 1, 0),
round_number = as.integer(round_number),
notResponded = ifelse(responded > 0, 0, 1) %>% as.integer()) %>%
arrange(round_number, ID, Group) %>%
drop_na()
)
This leaves us with 716 possible responses to requests.
We then fit the varying intercept and slope model, again with participants nested within groups.
b2.5 <-
brm(data = d.trim, family = bernoulli,
notResponded ~ 0 + intercept + round_number + Condition
+ (0 + intercept + round_number + Condition | Group/ID),
prior = c(prior(normal(0, 1), class = b),
prior(student_t(3, 0, 10), class = sd),
prior(lkj(1), class = cor)),
iter = 2500, warmup = 1000, chains = 4, cores = 4,
sample_prior = TRUE,
control = list(adapt_delta = 0.99),
seed = 2113)
save(b2.5, file = "models/exp2_brmsfit5.rda")
The priors we used.
prior_summary(b2.5)
The results.
print(b2.5)
## Family: bernoulli
## Links: mu = logit
## Formula: notResponded ~ 0 + intercept + round_number + Condition + (0 + intercept + round_number + Condition | Group/ID)
## Data: d.trim (Number of observations: 716)
## Samples: 4 chains, each with iter = 2500; warmup = 1000; thin = 1;
## total post-warmup samples = 6000
##
## Group-Level Effects:
## ~Group (Number of levels: 40)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept) 1.03 0.46 0.15 2.00 1.00 992 1023
## sd(round_number) 0.09 0.04 0.01 0.18 1.01 637 1085
## sd(Condition) 1.17 0.57 0.11 2.35 1.00 592 1054
## cor(intercept,round_number) 0.10 0.42 -0.67 0.87 1.01 814 1978
## cor(intercept,Condition) -0.39 0.41 -0.92 0.63 1.01 720 1881
## cor(round_number,Condition) 0.36 0.41 -0.63 0.93 1.01 703 1341
##
## ~Group:ID (Number of levels: 78)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept) 0.54 0.41 0.02 1.53 1.00 751 1393
## sd(round_number) 0.03 0.03 0.00 0.09 1.00 1445 2133
## sd(Condition) 0.58 0.46 0.02 1.68 1.00 708 1496
## cor(intercept,round_number) -0.09 0.50 -0.91 0.84 1.00 2882 3262
## cor(intercept,Condition) -0.44 0.50 -0.98 0.72 1.00 1071 2548
## cor(round_number,Condition) -0.05 0.51 -0.90 0.87 1.00 3727 4161
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept -1.70 0.36 -2.45 -1.02 1.00 1879 2954
## round_number -0.06 0.04 -0.15 -0.00 1.01 1360 2272
## Condition 0.49 0.39 -0.31 1.23 1.00 1919 2638
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Plotting the parameters.
stanplot(b2.5)
Let’s look at the trace plots to make sure Stan sampled efficiently.
plot(b2.5, ask = F)
Rhat values, n_eff values, and trace plots look okay.
The fixed effect of condition is 0.49, with 95% credible intervals that cross 0. This implies that participants were no more likely to not respond to requests in the Hidden condition.
As before, we sample from the posterior and convert this difference to the absolute probability scale.
post <- posterior_samples(b2.5)
visible_prob <- inv_logit_scaled(post$b_intercept)
visible_prob %>%
median() %>%
round(2)
## [1] 0.16
hidden_prob <- inv_logit_scaled(post$b_intercept + post$b_Condition)
hidden_prob %>%
median() %>%
round(2)
## [1] 0.23
difference <- hidden_prob - visible_prob
quantile(difference,c(0.025,0.5,0.975)) %>% round(2)
## 2.5% 50% 97.5%
## -0.05 0.08 0.19
The probability difference between the conditions is 0.07 (median), with 95% CIs crossing 0.
We compute a Bayes factor for this difference between probabilities.
(1 / hypothesis(b2.5, "inv_logit_scaled(intercept) - inv_logit_scaled(intercept + Condition) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 0.67
This Bayes factor implies anecdotal support for the hypothesis that the probabilities are equal across conditions.
bayes_R2(b2.5) %>% round(2)
## Estimate Est.Error Q2.5 Q97.5
## R2 0.24 0.03 0.18 0.3
bayes_R2() tells us that this model explains around 24% of the variance in the outcome.
Again, the data wrangling for this model is a little trickier.
(
d.trim <-
d %>%
# 1. select herd size, received, and requested
select(ID, Group, Counterbalancing, ends_with('.player.herd_size_after_shock'),
ends_with('.player.received'),
ends_with('.player.request_amount')) %>%
# 2a. flip the variables to get what the player gave and what their partner requested
mutate_at(vars(ends_with(".player.received"), ends_with(".player.request_amount")),
function(x) x[1:nrow(d) + c(1,-1)]) %>%
gather(key, value, -ID, -Group, -Counterbalancing) %>%
extract(key, c("Condition", "round_number", "variable"),
"SurvivalGame_(Hidden|Visible)[.](.|..)[.]player[.](herd_size_after_shock|received|request_amount)") %>%
spread(variable, value) %>%
# 2b. rename the variables to match their new meaning
rename(gave = received,
partner_requested = request_amount) %>%
mutate(Condition = ifelse(Condition == 'Hidden', 1, 0),
round_number = as.integer(round_number)) %>%
arrange(round_number, ID, Group) %>%
# 3. drop NAs, as no requesting happened
drop_na() %>%
# 4. was the player able to give?
filter(herd_size_after_shock - partner_requested >= 64) %>%
# 5. code 0 = request fulfilled, 1 = request not fulfilled
mutate(notFulfilled = ifelse(gave >= partner_requested, 0, 1) %>% as.integer())
)
This leaves us with 473 possible response decisions in which the player was able to give their partner what they asked for. Our outcome variable is whether they fulfilled that request or not.
We now fit the varying intercept and slope model, again with participants nested within groups.
b2.6 <-
brm(data = d.trim, family = bernoulli,
notFulfilled ~ 0 + intercept + round_number + Condition
+ (0 + intercept + round_number + Condition | Group/ID),
prior = c(prior(normal(0, 1), class = b),
prior(student_t(3, 0, 10), class = sd),
prior(lkj(1), class = cor)),
control = list(adapt_delta = 0.99),
sample_prior = TRUE,
iter = 2500, warmup = 1000, chains = 4, cores = 4,
seed = 2113)
save(b2.6, file = "models/exp2_brmsfit6.rda")
Here are the priors we used for model b2.6.
prior_summary(b2.6)
Let’s see the results.
print(b2.6)
## Family: bernoulli
## Links: mu = logit
## Formula: notFulfilled ~ 0 + intercept + round_number + Condition + (0 + intercept + round_number + Condition | Group/ID)
## Data: d.trim (Number of observations: 473)
## Samples: 4 chains, each with iter = 2500; warmup = 1000; thin = 1;
## total post-warmup samples = 6000
##
## Group-Level Effects:
## ~Group (Number of levels: 40)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept) 1.20 0.57 0.13 2.34 1.00 1028 920
## sd(round_number) 0.17 0.07 0.04 0.33 1.00 889 1603
## sd(Condition) 0.78 0.59 0.03 2.19 1.00 1314 2728
## cor(intercept,round_number) 0.22 0.39 -0.57 0.88 1.00 778 1217
## cor(intercept,Condition) -0.03 0.49 -0.87 0.87 1.00 3004 3312
## cor(round_number,Condition) 0.22 0.47 -0.78 0.92 1.00 2221 4171
##
## ~Group:ID (Number of levels: 73)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept) 0.65 0.48 0.02 1.79 1.01 1125 2106
## sd(round_number) 0.08 0.05 0.00 0.19 1.00 1072 1991
## sd(Condition) 0.77 0.59 0.03 2.20 1.00 1016 1971
## cor(intercept,round_number) -0.18 0.49 -0.91 0.82 1.00 1601 2741
## cor(intercept,Condition) -0.27 0.50 -0.95 0.79 1.01 1913 3600
## cor(round_number,Condition) 0.08 0.50 -0.85 0.90 1.00 2175 3743
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept -0.54 0.38 -1.34 0.20 1.00 2724 3783
## round_number -0.13 0.06 -0.26 -0.04 1.00 2098 2585
## Condition 0.95 0.39 0.15 1.71 1.00 3861 3596
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Plotting the parameters.
stanplot(b2.6)
Let’s look at the trace plots to make sure Stan sampled efficiently.
plot(b2.6, ask = F)
Stan sampled efficiently.
The fixed effect of condition is 0.95, with 95% credible intervals above 0. This implies that participants were more likely to not fulfill requests (when able) in the Hidden condition.
On the absolute probability scale.
post <- posterior_samples(b2.6)
visible_prob <- inv_logit_scaled(post$b_intercept)
visible_prob %>%
median() %>%
round(2)
## [1] 0.37
hidden_prob <- inv_logit_scaled(post$b_intercept + post$b_Condition)
hidden_prob %>%
median() %>%
round(2)
## [1] 0.61
difference <- hidden_prob - visible_prob
quantile(difference,c(0.025,0.5,0.975)) %>% round(2)
## 2.5% 50% 97.5%
## 0.04 0.23 0.40
The probability difference between the conditions is 0.23 (median), with 95% CIs above zero.
We compute a Bayes factor for this difference between probabilities.
(1 / hypothesis(b2.6, "inv_logit_scaled(intercept) - inv_logit_scaled(intercept + Condition) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 9.75
This Bayes factor implies moderate support for the hypothesis that the probabilities differ across conditions.
bayes_R2(b2.6) %>% round(2)
## Estimate Est.Error Q2.5 Q97.5
## R2 0.44 0.03 0.37 0.5
This model explains around 44% of the variance in the outcome.
As with b2.4b, we test for order effects by including an interaction.
b2.6b <-
brm(data = d.trim, family = bernoulli,
notFulfilled ~ 0 + intercept + round_number + Condition*Counterbalancing
+ (0 + intercept + round_number + Condition | Group/ID),
prior = c(prior(normal(0, 1), class = b),
prior(student_t(3, 0, 10), class = sd),
prior(lkj(1), class = cor)),
control = list(adapt_delta = 0.99),
sample_prior = TRUE,
iter = 2500, warmup = 1000, chains = 4, cores = 4,
seed = 2113)
save(b2.6b, file = "models/exp2_brmsfit6b.rda")
Here are the priors we used.
prior_summary(b2.6b)
Let’s see the results.
print(b2.6b)
## Family: bernoulli
## Links: mu = logit
## Formula: notFulfilled ~ 0 + intercept + round_number + Condition * Counterbalancing + (0 + intercept + round_number + Condition | Group/ID)
## Data: d.trim (Number of observations: 473)
## Samples: 4 chains, each with iter = 2500; warmup = 1000; thin = 1;
## total post-warmup samples = 6000
##
## Group-Level Effects:
## ~Group (Number of levels: 40)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept) 1.15 0.56 0.14 2.33 1.00 1407 1783
## sd(round_number) 0.17 0.07 0.04 0.33 1.00 845 1431
## sd(Condition) 0.84 0.62 0.04 2.28 1.00 1249 2772
## cor(intercept,round_number) 0.26 0.38 -0.50 0.91 1.00 834 1808
## cor(intercept,Condition) -0.08 0.49 -0.89 0.85 1.00 2254 3167
## cor(round_number,Condition) 0.23 0.46 -0.74 0.92 1.00 2665 3676
##
## ~Group:ID (Number of levels: 73)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept) 0.68 0.49 0.03 1.83 1.00 1380 2641
## sd(round_number) 0.08 0.05 0.00 0.20 1.00 965 2131
## sd(Condition) 0.84 0.61 0.03 2.24 1.00 1125 2050
## cor(intercept,round_number) -0.19 0.50 -0.92 0.84 1.00 1449 2622
## cor(intercept,Condition) -0.30 0.49 -0.95 0.77 1.00 1877 3222
## cor(round_number,Condition) 0.08 0.50 -0.86 0.90 1.00 2420 3671
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept -0.77 0.45 -1.67 0.10 1.00 3641 3946
## round_number -0.13 0.05 -0.26 -0.04 1.00 2117 2799
## Condition 0.64 0.52 -0.40 1.63 1.00 4099 3584
## Counterbalancing 0.44 0.54 -0.65 1.49 1.00 4149 4248
## Condition:Counterbalancing 0.48 0.60 -0.68 1.67 1.00 4403 4408
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
Plotting the parameters.
stanplot(b2.6b)
The interaction parameter has 95% CIs that cross zero, indicating no interaction effect.
Converting the fixed effects to the probability scale:
post <- posterior_samples(b2.6b)
visible_VisibleFirst_prob <- inv_logit_scaled(post$b_intercept)
visible_VisibleFirst_prob %>%
median() %>%
round(2)
## [1] 0.32
hidden_VisibleFirst_prob <- inv_logit_scaled(post$b_intercept + post$b_Condition)
hidden_VisibleFirst_prob %>%
median() %>%
round(2)
## [1] 0.48
difference_VisibleFirst <- hidden_VisibleFirst_prob - visible_VisibleFirst_prob
quantile(difference_VisibleFirst,c(0.025,0.5,0.975)) %>% round(2)
## 2.5% 50% 97.5%
## -0.08 0.15 0.37
When participants play the Visible game first, the absolute probability difference between the conditions is -0.03 (median), with 95% CIs crossing 0. Participants were no more likely to request from their partner in the Hidden condition.
We compute a Bayes factor for this difference between probabilities.
(1 / hypothesis(b2.6b, "inv_logit_scaled(intercept) - inv_logit_scaled(intercept + Condition) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 1.43
This Bayes factor implies moderate support for the hypothesis that the probabilities are equal in each condition.
What about for when participants play the Hidden game first?
visible_HiddenFirst_prob <- inv_logit_scaled(post$b_intercept + post$b_Counterbalancing)
visible_HiddenFirst_prob %>%
median() %>%
round(2)
## [1] 0.42
hidden_HiddenFirst_prob <- inv_logit_scaled(post$b_intercept + post$b_Counterbalancing +
post$b_Condition + post$`b_Condition:Counterbalancing`)
hidden_HiddenFirst_prob %>%
median() %>%
round(2)
## [1] 0.69
difference_HiddenFirst <- hidden_HiddenFirst_prob - visible_HiddenFirst_prob
quantile(difference_HiddenFirst,c(0.025,0.5,0.975)) %>% round(2)
## 2.5% 50% 97.5%
## 0.03 0.26 0.46
When participants play the Hidden game first, the absolute probability difference between the conditions is +0.20 (median), with 95% CIs above 0. Participants are more likely to request from their partner when beneath the threshold in the Hidden condition.
We compute a Bayes factor for this difference between probabilities.
(1 / hypothesis(b2.6b, "inv_logit_scaled(intercept + Counterbalancing) - inv_logit_scaled(intercept + Counterbalancing + Condition + Condition:Counterbalancing) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 7.07
This Bayes factor implies moderate support for the hypothesis that the probabilities differ across conditions.
bayes_R2(b2.6b) %>% round(2)
## Estimate Est.Error Q2.5 Q97.5
## R2 0.44 0.03 0.38 0.5
This model explains around 44% of the variance in the outcome.
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 16299)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_New Zealand.1252 LC_CTYPE=English_New Zealand.1252 LC_MONETARY=English_New Zealand.1252
## [4] LC_NUMERIC=C LC_TIME=English_New Zealand.1252
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] forcats_0.4.0 stringr_1.4.0 dplyr_0.8.3 purrr_0.3.2 readr_1.3.1 tidyr_1.0.0 tibble_2.1.3
## [8] tidyverse_1.2.1 ggplot2_3.2.1 cowplot_1.0.0 brms_2.10.0 Rcpp_1.0.2 bayesplot_1.7.0
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-140 matrixStats_0.55.0 xts_0.11-2 lubridate_1.7.4 threejs_0.3.1
## [6] httr_1.4.1 rstan_2.19.2 tools_3.6.1 backports_1.1.5 R6_2.4.0
## [11] DT_0.9 lazyeval_0.2.2 colorspace_1.4-1 withr_2.1.2 tidyselect_0.2.5
## [16] gridExtra_2.3 prettyunits_1.0.2 processx_3.4.1 Brobdingnag_1.2-6 compiler_3.6.1
## [21] rvest_0.3.4 cli_1.1.0 xml2_1.2.2 shinyjs_1.0 labeling_0.3
## [26] colourpicker_1.0 scales_1.0.0 dygraphs_1.1.1.6 ggridges_0.5.1 callr_3.3.2
## [31] digest_0.6.21 StanHeaders_2.19.0 rmarkdown_1.15 base64enc_0.1-3 pkgconfig_2.0.3
## [36] htmltools_0.3.6 readxl_1.3.1 htmlwidgets_1.3 rlang_0.4.0 rstudioapi_0.10
## [41] shiny_1.3.2 generics_0.0.2 zoo_1.8-6 jsonlite_1.6 crosstalk_1.0.0
## [46] gtools_3.8.1 inline_0.3.15 magrittr_1.5 loo_2.1.0 Matrix_1.2-17
## [51] munsell_0.5.0 abind_1.4-5 lifecycle_0.1.0 stringi_1.4.3 yaml_2.2.0
## [56] pkgbuild_1.0.5 plyr_1.8.4 parallel_3.6.1 promises_1.0.1 crayon_1.3.4
## [61] miniUI_0.1.1.1 lattice_0.20-38 haven_2.1.1 hms_0.5.1 zeallot_0.1.0
## [66] knitr_1.25 ps_1.3.0 pillar_1.4.2 igraph_1.2.4.1 markdown_1.1
## [71] shinystan_2.5.0 reshape2_1.4.3 stats4_3.6.1 rstantools_2.0.0 glue_1.3.1
## [76] evaluate_0.14 modelr_0.1.5 vctrs_0.2.0 httpuv_1.5.2 cellranger_1.1.0
## [81] gtable_0.3.0 assertthat_0.2.1 xfun_0.9 mime_0.7 xtable_1.8-4
## [86] broom_0.5.2 coda_0.19-3 later_0.8.0 rsconnect_0.8.15 shinythemes_1.1.2
## [91] ellipsis_0.3.0 bridgesampling_0.7-2